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I.  Background 

OSPAC  (Optimal  Synthesis  Program  for  Automatic  Control)  is  a  digital 

computer  program  written  in  Fortran  IV,  which  concerns  itself  with  the 

stationary  linear  quadratic  Gaussian  optimal  control  problem.  This  problem 

can  be  outlined  as  follows:   Consider  a  system  described  by 

x  (t)  -  A  x(t)  +  B  u(t)  +  v_  w(t) 

£  (t)  -  C  x(t) 

where 

A     is  an  n  x  n  plant  matrix 

x(t)  is  an  n  x  1  state  vector 

B     is  an  n  x  p  control  matrix 

u(t)  is  a  p  x  1  control  vector 

Y_     is  an  n  x  t  disturbance  matrix 

w(t)  is  a  t  x  1  disturbance  vector 

y_(t)  is  a  q  x  1  output  vector 

C     is  a  q  x  n  output  matrix 

Here,  w(t)  is  a  vector  of  linearly  uncorr elated,  zero  mean  white  noise 
signals  with  Gaussian  amplitude  probability  distribution  functions.  The 
elements  of  w(t)  are  assumed  to  be  sample  functions  from  n  random  processes 
which  are  each  ergodic  and  are  jointly  ergodic.  The  covariance  matrix  for 
w(t)  is 

E  [w(t)  wT(t  +  t)]  =   F  6(t) 
where  s(t)  is  the  unit  impulse  function. 

The  measured  quantities  or  sensor  signals  are 

z(t)  =  H  w(t)  +  v(t) 


where 

z(t)  is  a  u  x  1  measurement  vector 

H    is  a  u  x  n  measurement  matrix 

v(t)  is  a  u  x  1  measurement  noise  vector 
The  elements  of  v(t)  are  assumed  to  be  sample  functions  from  P  random 
processes  each  of  which  are  ergodic  and  jointly  ergodic.  The  covariance 
matrix  for  v(t)  is 

E  [v(t)  vT(t  +  t)]  =  G  6(t) 

The  system  is  assumed  to  be  completely  controllable  and  completely 
observable.   It  is  desired  to  find  the  control  function  u(t)  which  mini- 
mizes the  quadratic  scalar  index  of  performance 

J  =  lim  ±  f   [y_T(t)  Q  Z(t)   +   uT(t)  R  u(t)]  dt 
T-»    o 

where 

Q  is  a  q  x  q  symmetric  output  cost  weighting  matrix  and  at 
least  positive  semidef inite 

R  is  a  p  x  p  symmetric  control  cost  weighting  matrix  and 
positive  definite 

The  solution  to  the  linear  quadratic  Gaussian  control  problem  can  be  out- 
lined as  follows: 

a.)  The  optimization  problem  can,  by  the  called  Separation  Theorem, 
be  broken  up  into  two  separate  problems,  an  optimal  control  problem  and 
an  optimal  estimation  or  filtering  problem. 

b.)  The  optimal  estimation  or  filtering  problem  generates  an  optimal 
estimate,  xft)  of  the  state  x(t).  This  estimate  is  optimal  in  the  sense  that 

lim  i  f  ±(t  )  x(t)  dt 
T^o    o 


is  minimized,  where  5T(t)  is  the  estimation  error  defined  as 

x(t)  =  x(t)  -  x(t) 
The  optimal  estimator  (or  Kalman  filter)  has  the  form 

x(t)  =  A  x(t)  +  B  u(t)  +  K  [z(t)  -  Hx(t)] 

The  estimator  gains  are  given  by 

T  -1 
K  =  P  H  G 

where  P  is  the  error  co variance  matrix 

T 
E  [x(t)  x(t  +  t)]  =  P  5(t) 

P  is  the  positive  definite  solution  to  the  steady  -  state  filter  matrix 
Riccati  equation 

T  T  T     -1 

A  P  +  P  A     +Y_Fv_PH     GHP  =  0 

c.)  The  optimal  control  problem  generates  an  optimal  control  law 
u(t)  which  is  a  linear  function  of  the  estimated  state 

u_(t)  =  -  L  x(t) 
where  L  is  a  p  x  n  optimal  controller  gain  matrix.  The  gain  matrix  L  is 
identical  to  the  one  obtained  by  solving  the  optimal  control  problem  with 
no  system  disturbance,  exact  state  information,  and  the  index  of  perfor- 
mance given  by 

00 

J  =  J  [  ZT(t)  Qy(t)+  uT(t)  R  u(t)]  dt 

o 
the  controller  gain  matrix  L  is  given  by 

-1  T 
L  =  R   B   S 


where  S  is  the  positive  definite  solution  to  the  steady-state  control 
matrix  Riccati  equation 

T      T  -IT 

-SA-AS-CQC+SBR   BS=0 

It  can  be  shown  that  the  state  covariance  matrix 

E  [x(t)  xT(t  +  t)1  =  (P  +  M)  6(t) 
where  P  is  the  solution  to  the  filter  matrix  Riccati  equation  and  M  is  the 
positive  definite  solution  to 

(A  -  B  L)  M  +  M  (A  -  B  L)T  +  K  G  KT  =  0 

In  addition  to  the  solutions  outlined  above,  it  can  be  shown  that  the 
transfer  matrix  relating  the  Laplace  transform  of  the  optimal  control  law 
u(t)  to  the  Laplace  transform  of  the  measurement  vector  z_(t)  (with 
v  (t)  =  C)  is  given  by 

U  (S)  =  -  L  (SI  -A  +BL  +  K  H)"1  K  Z  (S) 
where 

u(s)  -a [u(t)i 

Z  (S)  =•  J[z(t)l 

In  addition,  the  characteristic  roots  of  the  estimator  are  the  roots  of 

|SI-(A-KH)|=0 
and  the  characteristic  roots  of  the  state -feedback  controller  are  the  roots 
of 

|  S  I  -  (A  -  B  L)  |  =  0 

The  characteristic  roots  of  the  entire  closed-loop  system,  i.e.,  the  plant, 
estimator  and  state -feedback  controller  are  just  the  estimator  roots  and 
state  feedback  controller  roots  taken  together. 
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II.  OSPAC  Description 


A.  Introduction 

OSPAC  makes  extensive  use  of  the  Variable  Dimension  Automatic  Synthesis 
Program  (VASP)  configured  by  John  S.  White  and  Homer  Q.  Lee  of  NASA  Ames 
Research  Center.  A  documentation  report  entitled  "Users  Manual  for  the 
Variable  Dimension  Automatic  Synthesis  Program  (VASP),"  Oct.  1971,  may  be 
obtained  from  NTIS  (N72-10190).   OSPAC  can  provide  the  following  output: 

1.)  S,  the  solution  to  the  steady-state  control  matrix  Riccati 
equation. 


2. 
3. 

k. 

5. 
6. 

7. 
8. 
9. 


L,  the  controller  gain  matrix. 

P,  the  solution  to  the  steady-state  filter  matrix  Riccati  equation, 

K,  the  estimator  (filter)  gain  matrix 

P  +  M,  the  covariance  of  the  system  state 

CVHX,  the  covariance  of  H  x(t) 

CVU,  the  covariance  of  u(t) 

J,  the  value  of  the  index  of  performance 

the  roots  of 


|SI-(A-BL)|=0 

I  S  I  -  (A  -  K  H)  I  =  0 

10.)  the  elements  of  the  transfer  matrix  relating  U(S)  to_z(S) 
OSPAC  solves  the  steady-state  Riccati  equations  by  integrating  the 
differential  Riccati  equations  until  a  steady-state  solution  is  reached  or 
when  the  maximum  number  of  integration  steps  (as  specified  by  the  user,  see 
NC0NT(2)  below)  has  been  reached. 


B.  OSPAC  Input 

As  presently  configured,  the  maximum  dimensions  of  the  input  matrices 

for  OSPAC  are 

n  =  10 

p  =  10 

q  =  10 

t  =  9 

u  =  9 

It  is  imperative  that,  in  his  program,  the  user  ensure 

q  <  n 
t  <  n 
u  <  n 

If  the  conditions  above  are  not  met,  subroutine  AUG  will  produce  incorrect 

results. 

The  input  card  arrangement  is  shown  in  tabular  form  on  the  next  page. 

The  description  of  the  items  in  the  table  follows. 

NSOL  This  single  integer ,  in  format  (11)  specifies  the  number  of  problems 
to  be  run. 

NOPT  This  single  integer,  in  format  (il)  specifies  the  solution  option.   If 

NOPT  =  1,  one  obtains  only  the  state-feedback  controller  solution  (L,  S) 

NOPT  =  2,  one  obtains  only  the  estimator  or  filter  solution  (K,  P ) . 

NOPT  =  3>  one  obtains  the  controller,  estimator  and  covariance  solutions 
(L,  S ,  K,  P,  P  +  M,  CVHX,  CVU,  j) . 

NOPT  =  k,   one  obtains  the  same  solutions  as  in  NOPT  =  3  plus  the  system 
characteristic  roots  and  transfer  matrix. 
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NCONT  This  vector  of  length  3  is  input  for  each  Riccati  solution  (S,  P 
and  M)  in  format  (3110) 

NCONT(l)  =  1 

NCONT (2)  =  maximum  number  of  integration  steps  in  Riccati  solution; 
NC0NT(2)   should  be  ;>  100. 

NC0NT(3)  =  1 

A,  _B,  etc.   With  the  exception  of  NCONT,  each  input  matrix  requires  a 

header  card  in  format  (AU,  i+X,  2lU) .  This  represents  a  k 

character  title,  k   blank  spaces,  then  the  number  of  rows  and 

columns  in  the  matrix.  Each  row,  beginning  on  a  new  card  is 

entered  after  the  header  card.  Since  the  program  can  handle 

some  10  X  10  matrices  and  since  ihe  input  format  is  (7E10.5)j 

some  matrices  may  require  2  cards  per  row.  However,  each 

matrix  row  must  begin  on  a  new  card. 

C.  OSPAC  Output 

The  following  problems  may  occur  in  some  OSPAC  executions. 

1.)  UNDERFLOW  Messages;  the  VASP  programs  used  in  OSPAC  frequently 

generate  very  small  numbers  which  result  in  UNDERFLOW  error  messages.  The 

main  program  includes  an  ERRSET  subroutine  which  prevents  the  messages  from 

being  printed  each  time  such  an  "error"  occurs.  These  underflows  do  not 

compromise  the  solution  in  any  way. 

2.)  OVERFLOW  Messages;  If  OSPAC  generates  OVERFLOW  error  messages 

which  are  not  attributable  to  user  input  errors,  then  input  scaling  is 

necessary.  This  is  discussed  in  Section  III 
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3.)  Failure  to  converge;  Each  Riccati  solution,  S,  P  and  M,  is 
preceded  by  a  statement  indicating  the  number  of  iterations  required  to 
obtain  the  solution.   If  this  number  equals  the  value  of  NC0NT(2)  for  that 
equation,  then  the  solution  has  not  converged.  Assuming  that  the  system 
is  controllable  and  observable,  and  that  NC0NT(2)  ^  100,  failure  to  con- 
verge usually  means  that  the  integration  step  size  is  too  large.  The  VASP 
routines  are  supposed  to  automatically  adjust  the  initial  stepsize  (set 
equal  to  1.0D+00  in  the  third  argument  of  the  subroutines  ETPHI  called 
in  the  main  program)  for  each  problem.  Occasionally,  however,  this  auto- 
matic procedure  fails.  The  user  should  then  reduce  the  value  of  the 
constant  in  the  ETPHI  call  statement  for  the  particular  equation  (controller, 
filter,  or  state  covariance)  which  failed  to  converge,  and  rerun  the  job. 

The  appendix  provides  a  listing  of  the  main  program  and  all  the 
subroutines. 

III.  Scaling  Considerations 
For  large  systems  (n  ^  8),  some  scaling  of  the  input  matrices  is 
often  necessary.  When  OSPAC  produces  overflow  error  messages  and  no  ob- 
vious source  for  these  errors  can  be  found,  scaling  is  indicated.  A 
simple  scaling  procedure  that  has  been  used  with  a  good  deal  of  success 
with  OSPAC  involves  amplitude  scaling  the  system  equations  as  though 
they  were  going  to  be  programmed  on  an  analog  computer.  Again  consider 
the  system 

x  (t)  =  A  x(t)  +  B  u(t)  +  ^  w(t) 

jr  (t)  =  C  x(t) 

z  (t)  =  H  x(t)  +  v(t) 

m 

J  =  11m  i  J*  (x(t)  Q  v_(t)  +  uT(t)  R  u(t)]  dt 
T_  a    o 
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Now  define  the  following  matrices: 

Xw  is  an  n  x  n  diagonal  matrix  whose  elements,  x   ,  are  "guestimates" 
— M  m .  . 

of  the  maximum  values  of  the  state  variables  x.(t). 

l  ' 

U„  'is  a  p  xp  diagonal  matrix  whose  elements,  u._  ,  are  ''guest imates"  of 
~"M  M.  . 

11 

the  the  maximum  values  of  the  controls  u.(t). 

1  ' 

Yw  is  a  q  x  q  diagonal  matrix  whose  elements,  yw  ,  are  "guestimates"  of 
— M  M .  . 

n 

the  maximum  values  of  the  output  variables  y.(t). 

Zw  is  a  u  x  u  diagonal  matrix  whose  elements,  z„  ,  are  "guestimates"  of 
— M  M.  . 

n 

the  maximum  values  of  the  measurements  z.(t). 

l 

Now  define 

£s(t)     =     *£  *(*) 

£s(t)    -     tftft) 

ajt*)   ■   £**■> 

where  the  subscript    '  s'   refers  to  scaled  quantities.     Rewriting  the  original 
system  equations  using  the  matrices  defined  above  yields 

W*)   -£^xs(t) 

Z^t)      -HX^t)    +    v(t) 

J="'f    I    [Cym  ys(t)^    4  [iH  Jr,(t)l  +CUM  «s(t)]T  R 

ir*  oo  o  —  —  ■*— 

[UMug(t)]}  dt 
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These  equations  can  be  written 


x  (t)  =  A  x  (t)  +  B  u  (t)  +  y  w(t) 

sv  '     s  sv  '  s  s       s  — N  ' 


ys(t)  =CsXs(t) 


z  (t)  =  H  x  (t)  +  u  (t) 
sv  '  s  sx  '  s 


with 


where 


T 
J  =  lim  i  f  [y  T(t)  0  y  (t)  +  u  T  (t)  R  u  (t)]  dt 


F  =  F 
s   — 

G  =  Z"1  G  Z,,"1 

s    M  —  M 


-1 

"M 


A-  =  Xjj  AX 


B  =  X'1  B  II, 
s    M  —  M 


s    M  —  M 


Ys  ■  ^  1 


H  —  Z.,  H  X,. 

s    M  —  M 


Q  =  Y„T  Q  YM 

s    M  —  M 


R   =  U  J   R  U 

s    M  —  M 
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The  scaled  matrices  above  are  then  used  as  inputs  to  OSPAC.  The 
output  of  OSPAC  can  then  be  unsealed  to  obtain  the  solution  to  original 
problem.  Unsealing  the  pertinent  output  quanties  is  summarized  below. 

T 

p  =  Xm  p  x™ 

—    M  s  M 


L  =  LXM 


-1 


P  +  M  =  XM  (P  +  M)  XMT 

—   —   _M  —   —  s  _M 

SB)   =^(S) 

=  U  [-L   (SI  -  A  +  B  L  +  K  H  )_1  K  ]Z  "1  2(S) 
M   — s   —   _s   _s_s   _s  _s    _s.  _M 

It  should  be  emphasized  that  the  eigenvalues  of  the  problem  are  unaffected 
by  amplitude  scaling,  i.e.  the  roots  of 

|  S  I  -  (A  -  B  L)  |  =  0 

and 

I  "  I  "  (*,  -  B,  L,)  |  =  0 

are  identical,  as  are  the  roots  of 

|  S  I  -  (A  -  K  H)  |=0 
and 

|  S  I  -  (A  -K  H  )  I  =  0 
1   -   s-s   -s  —s.l 

as  are  the  roots  of 

|SI-A+BL  +  KH|=0 


and  |SI-A+BL+KH|=0 

—       — s  s  — s        — s  -s    I 


Ik 


IV.  Sample  Problem  -  Helicopter  Optimal  Control  Problem 

The  longitudinal  motion  of  a  helicopter  near  hover  in  turbulence  can 

be  modeled  reasonably  well  by  the  following  set  of  differential  equations 

••     • 

9  -  a-|Q  -  a2  (u  -  u  )  -  b6  =  0 

•     • 

u  -  a  6  -  a^  (u  -  u  )  -  g(e+6)  =  0 


where 


u  =  longitudinal  fore -aft  turbulence  (here  assumed 
to  have  a  white  spectrum),  ft/sec 

9  =  pitch  angle  of  fuselage,  rad 

6  =  tilt  of  rotor  tip  path  plane  with  respect  to  fuselage, 
rad 

o 
g  =  acceleration  due  to  gravity,  ft/sec 

u  =  groundspeed  measured  from  trim,  ft/sec 


vertical 


rotor  thrust  line 


fuselage  reference  axis 


-^u 


The  synthesis  problem  centers  about  finding  the  control  law  6(t)  which 
minimizes  a  quadratic  index  of  performance  with  groundspeed  being  the 
measured  variable. 
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a     =  -   .k  /sec 


a     =  -  U.593  ft/sec 


a     =  -    .0030l+8/ft-sec 


a.    =  -     .02/sec 


Assume 


E[u   (t)  u  (t  +  t)] 


b       =    -    6.3/sec' 


25  6(T)   ft2/sec2 


The  measured  variable  is     u     and 

E[v(t)  v(t  +  t)]    = 
One  can  define  a  set  of  state  variables 


,01  s(T)  ft2/sec2 


xl 

e 

* 

X2 

— 

e 

X3 

= 

u 

and  a  set  of  state  equations 


A 


Now 


L 


0 
0 

g 


►  = 


xl 

F 

2 

1X2 

►         + 

b 

if 

_ 

^ 

_g. 

6  + 


z  =         [0,  0,   lj 


,    0,    11   <  X.  >  +      V 


1  0 
0  1 
0         0 


LX3J 


0 
0 

1 


w 


v.    3' 


-a. 


-a 


LXJ 


u 


g 
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Q  = 


100. 
0 
0 


0 

100. 

0 


0 
0 
.0U 


R  =  100. 


Note  that,  if,  at  some  instant  of  time 


9  = 

e  = 


u  = 


6  = 


.1  rad 

.1  rad/sec 
5.0  ft/sec 
.1  rad 


each  scalar  term  in  the  integrand  of  the  index  of  performance  would  he 
making  a  contribution  of  unity  to  the  integrand.  This  provides  some 
rationale  for  the  selection  of  the  Q,  and  R  matrices  above. 

The  input  deck  set-up  is  shown  on  the  next  pages.  Following  that  is 
the  0SPAC  output.  No  scaling  was  necessary  in  this  problem. 
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Helicopter  Optimal  Control  Problem 
Input  Deck 


1 

?  |   5  |   «  |   » 

6 

7   |    8   |    9  |  10 

column 

ii|iJM3|i41i5|i6|i7|ie|i9|ZOI?i|?2|23|2«|i5j26|a7|28|29|30IJI|J2|S3|34|J5|J6|37|je!39|40 

•' 

/ 

1        1        1       1       1        1       1       1       1 

1     1     1     1     1     i     !     1     1 

H 

EiLi  iiC 

rf 

Pi71E:K 

i(*;Pmiirti>kiLi    \C 

0.AJlT.Cifl4lLi     iPiKid  BiL  llcirti     l     i     i     I     i 

4 

1    1    1    1   1    1    1    1    1 

1     1     1     1     l     I     1     1     1 

/ 

/1O1O1    i    i    i    i    i    i 

/,,,,.,, 

A 

i3i    i    i    i3i    i    i    i 

I     1     I     I     l     i     I     1     i 

0. 

/ 1- 1    i    i    i    i    i    i    i 

01.  1     1     1     1     I     I     I     1 

0 

•  i    i    i 

""  i  - 1 ^i    i    i    i    i    i    i 

-i.  i6i0i3!0i<?i8i    i 

3 

2i.  i2i 

-i-?i-i$-|7i3i    i    i    i 

-i.i<)i.2i    i    i    i    i    i 

8 

i3i    i    i    i/i    i    i    i 

i    i    i    i    i    i    i    i    i 

0 

.  i     ii 

i    i    i    i    i    i    i    l    i 

1       1      1      1      1      !      1      1      1 

<e 

.  i3i    , 

i    i    i    i    i    i    i    i    i 

1      1      1      1      1      1      1      1      I 

3 

-?  i.  r-ii 

i    i    i    i    i    t    i    i    i 

1      I      1      1      1      I      1      1      1 

C 

i3i    i    i    i3i    i    i    i 

1      1      1      1      1      1      1      1      1 

/ 

*  1      1      1 

Oi . i    i    i    i    i    i    i    i 

Oi..     i     i     i     i     i     i     | 

- 

0 

*  i       1 

i    i    i 

l\  .i    i    i    i    i    i    i    i 

Oi  .  i    i     i    i    i    i    i     i 

6 

Oi.i    i    i    i    i    i    i    i 

/i  .1    i    i    i    I    i    i    I 

<k 

i3i    i    i    1 3i    i    i    | 

i     i    i     i    i    i    I    i    i 

1 

OlOl.i 

Oi .  i    i    i    i    i    i    i    i 

Oi.i    i    i    i    i    i    i    i 

b 

.  i     i     i 

/   I  0\  Ol  .  1       1       1       1       1       1 

tfi.i  i  i  i  i  i  i  i 

0 

«i     i     i 

0\.  i    i    i    i    i    i    i    i 

.  idi^i  i 

i 

2  i J  i  *  i i 

6 

y  i  •  |  9  |  io 

1 1  |  12  j  13(14  |  IS  1 16  |  17  |  II  119120 

21  |22| ?3|24|25| 26i27|2S(?9j  30 

34  132|33|Mi33i36|37l3«|3t|40 

4i 
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Input  Deck   (cont'd) 


1 

3 1 s  i  *  i 5 

E 

7    j    8   |    9   |  iO 

column 

ii|:;liJ|i4liJ|i6|i7|i«|i9|?0,Ji|?2l25l2«|^il26|27l28|29|JOlSllJ2|J3|J4|35|36lir|J«!39|40|«i 

R 

i    i    i 

,/,        .../,.,. 

1 

0  O  i .  i 

1 

1    1    1 

/i0i6i    i   i    i    i    i    i 

/,                  1         1         1         1         1        I         1 

F 

1    1    1 

i/i    i    i    i/i    i    i    i 

a. 

51.1    i 

i    '  i       i              it                      i 

6 

i    i    i 

i/i    i    i    i/i    i    i    i 

• 

0./,    i 

0 

A,H,    , 

i3i    i    .    i/i    i    .    i 

c 

-iii 

I       1       1        1        1       1       I       1        1 

* 

0,6,3,0 

1 

&     I     , 

r 

Oi^i    i 

I       I        i       I        I        i       I        1       1 

1      I      1      1      1      !      I      1     1 

H 

i    i    i 

i/i    i    i    i3i   ■    i   i 

0 

•  i    i    i 

Oi.i    i    i    i    i    i    i    i 

/■-. 

1 

l    i    i 

/lOlOl      1      1      I      1      1      ! 

/,      ,      .      .      .      1      1      .     1 

i    i    i 

i    i    i 

i    i    i 

i    i    i 

i    i    i 

i    i    i 

, 

t|)|4|S 

6 

7   ]    •    |   9  |  iO 

II  |  12]  13}  14  |  19  J  16  |  i  7  |  18  |  I9|?0 

21  |22|23|24|25|26|2T|2»|29|30 

34    |32|33|34|  3S|36(  JT|3«|3»|40 

41 
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Output 


HELICOPTER    OPTIMAL    CONTROL    PROBLEM 


A      MATRIX  3    ROWS  3    COLUMNS 

0.0  l.OOOOOOOD    00         0.0 

0.0  -4.0000000D-01  -3.04800000-03 

3.22000000    01       -4.59300000    00  -2.00000000-02 

B       MATRIX  3    ROWS  1    COLUMNS 

0.0 

6.3000000D  00 
3.22000000    01 

C      MATRIX  3    ROWS  3    COLUMNS 

l.OOOOOOCD    00  0.0  0.0 

0.0  1.00000000    00  0.0 

0.0  0.0  l.OOOOOOOD    00 

Q       MATRIX  3    ROWS  3    COLUMNS 

l.OOOOOOOD    02         0.0  0.0 

0.0  1.00000000    02  0.0 

0.0  0.0  4.00000000-02 

R      MATRIX  1     ROWS  1    COLUMNS 

l.OOOOOOOD    02 

10  ITERATIONS 

S               MATRIX                    3    ROWS  3    COLUMNS 

1,94617480    02          1.2902733D    01  2.48101S5D    00 

1.2902733D    01          1.79268850    01  -  2.0 154349D-01 

2.4810185D    00      -2 .0 1 54349D-0 1  9.93  738  860-02 

L  MATRIX  1    ROWS  3    COLUMNS 

1.6117602D    00         1.06449680    00         1 .9301 1 5 1D-02 

F      MATRIX  1    ROWS  1    COLUMNS 

2.5OOOO00D    01 

G      MATRIX  1     ROWS  1    COLUMNS 

1.0000000D-02 

GAM      MATRIX  3    ROWS  1    COLUMNS 

0.0 

3.0480000D-03 
2.0000000D-02 

H      MATRIX  1    ROWS  3    COLUMNS 

0.0  0.0  l.OOOOOOOD    00 

12         ITERATIONS 

P  MATRIX  3    ROWS  3    COLUMNS 

9.1510898D-05  7 .84982 78D-05  1. 25297940-03 

7.8498278D-05  1 . 59 593240-04  9.9  2632430-04 

1. 25297940-03  9 . 9263243D-04  2.83617550-02 

K  MATRIX  3    ROWS  1    COLUMNS 

1.2529794D-01 
9.92632430-02 
2.8361755D    00 

10  ITERATIONS 

P  +  M         MATRIX                     3    ROWS  3    COLUMNS 

2.0091399D-04         4 . 38083 5  ID- 10  -4.5  0784550-04 

4.38083510-10         3 . 7590S46D-04  -4.55009350-03 

4.5078455D-04      -4. 55009 35D-03  4. 79222 65D-01 
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Output   (cont'd) 


CVU         MATRIX  1     POWS  i    COLUMNS 

9.4159189D-05 

CVHX       MATRIX  1     ROWS  1    COLUMNS 

4.7922265D-01 

THE    INOEX    OF    PERFORMANCE,     J:  0.086 

THE    ZEROS    OF:    DET (S I- ( A-BL ) ) 

REAL         IMAGINARY 

-0.62793D  01     0.0 
-0.734280  00    -0.32735D  00 
-0.734280  00     0. 32735D  00 

THE  ZEROS  OF:  DET ( S I- ( A-KH) ) 

REAL         IMAGINARY 

-0.21281D  01     0.0 
-0.56406D  00    -0.I4101D  01 
-0.564060  00     0.14101D  01 

THE  ZEROS  OF:  THE  DENOMINATOR  POLYNOMIAL  OF  U(S)/Z(S) 

REAL         IMAGINARY 

-0.86889D  01     0.0 
-0.94757D  00    -0.25163D  01 
-0.94757D  00     0.25163D  01 

THE  COEFFICIENTS  OF  THE  POLYNOMIAL  IN  INCREASING  POWERS  OF  S 

0.62816D  02    0.23696D  02    0.10584D  02    0.10030D  01 

THE  ZEROS  OF:  THE  NUMERATOR  POLYNOMIALS  OF  U(S)/Z(S) 

1=   1      J=   1 

THE  ZEROS  OF  THIS  ELEMENT 

REAL         IMAGINARY 

-0.22480D  01    -0.39099D  01 
-0.22480D  01     0.39099D  01 

THE  COEFFICIENTS  OF  THE  POLYNOMIAL  IN  INCREASING  POWERS  OF  S 

0.11528D  04    0.25481D  03    0.56675D  02 
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APPENDIX 


Program  Listing 
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OSPAC 
OPTIMAL  SYNTHESIS  PROGRAM  FOR  AUTOMATIC  CONTROL 

THIS  PROGRAM  MAKES  EXTENSIVE  JSE   OF  THE  VARIABLE   DIMENSION 
AUTOMATIC  SYNTHESIS  PROGRAM  (VASP)  CONFIGURED  BY  JOHN  S.   WHITEt  ET 
AL.t  AMES  RESEAKCH  CENTER,  OCTOBER  1971.  A  DOCUMENTATION  REPOPT 
(N72-10190)  MAY  3E  03TAINED  FROM  THE  NT  I S  FOR  VASP. 

OSPAC  SOLVES  THE  FILTER,  CONTROL,  AND  STATE  COVARIANCE  MATRIX 
RICCATI  EQUATIONS  FOR  A  SYSTEM  DEFINED  BY: 

X=AX+BU+GAMW 


Z=HX*V 

IN  ADDITION,  THE  FILTER  AMD  CONTROL  CHARACTERISTIC  ROOTS  ARE  FOUND 
ALONG  WITH  THE  TRANSFER  MATRIX  U(S)/Z(S) 

INPUT  MATRICES 

A  N  BY  N  PLANT 

B  N  BY  P  CONTROL 

C  Q  BY  N  OUTPUT 

F  T  BY  T  DISTURBANCE  COVARIANCE 

G  U  BY  U  MEASUREMENT  NOISE  COVARIANCE 

H  U  BY  N  OBSERVATION 

R  P  BY  P  CONTROL  COST  WEIGHTING 

Q  Q  BY  Q  OUTPUT  COST  WEIGHTING 

GAM  N  BY  T  DISTURBANCE 

OUTPUT  MATRICES 

P  N  BY  N  ERROR  COVARIANCE   (FILTER  Ri:CATI  SOLUTION) 

S  N  BY  N  (CONTROL  RICCATI  SOLUTION) 

M  N  BY  N  (COVARIANCE  INTERMEDIATE  SOLUTION) 

P+M  N  BY  N  STATE  COVARIANCE 

K  N  BY  U  OPTIMAL  ESTIMATOR  GAIN  (FILTER) 

L  P  BY  N  OPTIMAL  CONTROLLER  GAIN  (CONTROL) 

CVHX  U  BY  U  COVARIANCE  OF  HX 

CVU  P  BY  P  COVARIANCE  OF  U 

SYSTEM  VECTORS 

X  N  BY  1  STATE 

U  P  BY  1  INPUT 

Y  Q  BY  1  OUTRUT 

V  U  BY  1  MEASUREMENT  NOISE 

W   T  BY  1    SYSTEM  DISTURBANCE 
Z   U  BY  1    MEASUREMENT 

OUTPUT  SCALARS 

J   INDEX  OF  PERFORMANCE 

SYSTEM  CHARACTERISTIC  ROOTS  AND  ELEMENTS  OF  U(S)/Z(S)  TRANSFER 
MATRIX 

SCALAR  SOLUTION-CONTROL  PARAMETERS 

NSOL   NUMBER  OF  ^ROBLE^S  TO  BE  RUN   (II) 
NOPT   TYPE    OF  PROBLEM  TO  BE  RUN    (II) 

=1   CONTROL  ONLY  (L,S) 

=2   FILTER  ONLY  (K,P) 

«=3   EVERYTHING    (  L  ,  S  ,  K  ,  P,  M  ,  MP,  Z  1  ,  J  ) 

=4   SAME  AS  3  BUT  WITH  CHARACTERISTIC  ROOTS  AND  TRANSFER 
MATRIX 


A  1 


C  VECTOR    SOLUTION-CONTROL    PARAMETERS 

C 

C  NCONT(I),  1=1,2,3   (3110) 

C  NC0NT(1)=1 

C  NC0NT(2)=THE  MAXIMUM  NUMBER  OF  STEPS 

C  NC0NT(3)=1 

C 

C  WITH    THE    EXCEPTION!    OF    NCONT,     EACH    MATRIX    INPUT    REQUIRES    A    HFADFR 

C       CARD    OF    FORMAT     (A4,4X,2I4)     (A    ^-CHARACTER    TITLE,    4    BLANKS,     AND    THE 

C       NUMBER    OF    ROWS    AND    COLJMNS    IN    THE    MATRIX).       EACH    ROW,     BEGINNING    ON 

C       EACH    ROW.       HOWEVER,     EACH    MATRIX    ROW    MUST    BEGIN    ON    A    NEW    CARD. 
c __ 

C       THE    FORMAT    FOR    INPUT    MATRICES     IS     (7E10.5) 

C 

C 

C  MAXIMUM    DIMENSIONS: 

C 

C  N=10           NOTE:   THE  USER  MUST  ENSURE  THAT  Q,T,  AND  U 

C  P=10                   ARE  EACH  LESS  THAN  N  IN  HIS  PROGRAM 

C  Q=9 

C  T=9 

C  U=9 

C 

C  INPUT  CARD  DECK  ARRANGEMENT 

C 

C        CARD  # 

C  1    NSOL   (U) 

C  2    TITLE  FOR  PROBLEM  (A72) 

C  3    NOPT   (II) 

C  4+ 

C  MATRICES  REQUIRED  <ZZ  CARDS) 

C  N0PT=1       NCONT,A,B,C,Q,R 

C  N0PT=2       NCCNT,A,F,3,GAM,H 

C  N0PT=3       NCONT , A, B,C,Q,R, NCONT, F,G, GAM, H, NCONT 

C  NOPT  =  <*       SAME    AS    3 

C  4  +  ZZ     TITLE  FOR  NEXT  PROBLEM 

C  ETC. 

C 

C 

C         

C 
C 
C 
C 

c 

Q  *************** 

c  *  * 

C  *   DIMENSION   * 

C  *  * 

C  *************** 

c 

DOUBLE  PRECISION   A(  10  , 1 0 ) , AT( 1 0  , 1 0 }  ,  B ( 10 , 10 ) , C ( 10, 1 0 ) , D ( 1 0,  10 ) , 

1  3<10,10),H(10,10),R(10,13),Q(10,10) ,GAM( 10,10) , 

2  P(10,10),S(10,10),M(10,10),MP(10,10),Z1(10,10), 

3  <(  10,10)  ,L(  13,10)  ,Z (23,23)  , PHI  (23,20  )  , 

4  PHIT(20,20),F(13,10),RI(13,10),GI(10,10), 

5  DUM(100  0),CT(10,10),J,TRA,TRB,GAMT(10,10), 

6  HTUO, 10), KT(10, 10),  J  (23), V(23),W3(10),W4  (10), 

7  ZED(10,10,10),CH(20),DC(10,10),LT(10,10), 

8  MLTdOt  1J)  tCVU(lOilO) 
C 

DIMENSION  NA(2) , NAT ( 2 ) , NB ( 2 ) ,NC(2) ,NO(2),NG(2) ,NH(2) , NR ( 2 ) ,NQ(2 ) , 

1  NGAMC2) ,MP(2) ,NS(2) , NM(2)  ,NMP(2)  , NZ 1 ( 2 ) , NK ( 2 ) ,ML ( 2  ) , 

2  NZ(2),^DHI(2),MPHIT(2),NFi2),NRI(2),NGI(2),NCT(2), 

3  NCONT(3),NGAMT(2),NHT(2) ,NKT(2) , NLT ( 2 ) , NML T ( 2 ) , NC VU( 2 ) 
C 

CALL     ERRSET     (203,300,-1,1) 
COMMON       /MAX/MAXRC 
COMMON    ZED 

C0MM0N/LINES/NLP,LIM,TITLE(23) 
KDUM    =    1000 


A   2 


MAXRC  =  100 
C   HOW  MANY  SOLUTIONS? 

READ     (5tl3)     NSOL 
C 

DO    12    NTIMES=1,NS3L 

CALL    RDTITL 

READ     (5,13)     NOPf 
C 

c 

c  *  * 

C  *      CONTROLLER    SOLUTION      * 

C  *  * 

c  *$*$*$*************#***** 

c 
c 

IF    (N0PT.EQ.2)    GO    TO    I 

READ    (5,14)     (NCONlTt  I)  ,1=1  ,3) 

CALL    READ    ( 5 , A ,NA, B , NB,C , NC,Q , NQ, R , NR ) 

NS(1)     =    NA( 1) 

NS(2)  =  NA( 2) 
C  FIND  R  INVERSE 

CALL  EQUATE  ( R,NR, RI , NRI ) 

CALL  INV  (RI,NIRI,DET,DUM) 
C  SET  S=0 

CALL  SCALE  ( S , NA, S , NS , 0. DOO ) 
C  FIND  EXPONENT,  FIND  RI*3  TRANSPOSE 

CALL  AUG  (A ,NA,6,N3,RI ,NRI , C , NC , Q, NQ , D, ND, Z , NZ , 0 ) 
C  FIND  EXP  Z 

CALL  ETPHI  (Z,NZ,1 . ODOO, PH I ,NPH I ,DUM, KDJM ) 
C  FIND  RICCATI  SOLUTION 

CALL  RICAT  ( PHI , NPHI ,D ,ND ,NCONT ,L , NL , S,NS , DUM, KDUM) 

CALL  PRNT  (S,N5,4HS    ,1) 

CALL  PRNT  (L,NL,4HL    ,1) 
C  CHECK  NOPT  TO  SEE  IF  3D\JE 

IF  (NOPT.EQ.l)  GO  TO  12 
C  ELSE,  DO  FILTER  SOLUTION 
C    READ  ADDITIONAL  CONTROL  AND  MATRIX  CARDS  REQJIRED 

READ  (5,14)  (NCONT( I ) ,1=1 ,3) 

CALL  READ  (  4,F  ,MF  ,  3  ,  NG  ,3AM  ,N3AM  ,  H,  N,H  ,H  ,NH) 

GO  TO  2 
C 

c  *  * 

C  *   FILTER  SOLUTION   * 

C  *  * 

C 

c 

1  READ  (5,14)  (NCONT(  I  ),  1  =  1  ,3) 

CALL  READ  ( 5 , A , NA, F , NF , G, NG,GAM ,H, NH ) 
C  FIND  A  TRANSPOSE 

2  CALL  TRAMP  ( A ,NA , AT , NAT ) 

CALL  TRANP  ( G AM, NGAM, GAMT , NGAMT ) 

CALL  TRANP  ( H , NH , HT , NHT ) 
C  FIND  G  INVERSE 

CALL  EQUATE  ( G  ,N3 ,GI , NGI ) 

CALL  INV  (GI,NGI ,DET,DUM) 
C  SET  P=0 

CALL  SCALE  ( P , NA , P, MP , 0. DOO) 
C  FIND  EXPONENT,  G  INVERSE*  H 

CALL  AUG  ( AT, NAT, HT, NHT, GI, NGI, 3AMT, NGAMT, F,NF,D,ND,Z, NZ,0) 
C  FIND  EXP  Z 

CALL  ETPHI  (Z, NZ , 1 . DOO, PHI , NPHI , DUM, KDUM) 
C  FIND  RICCATI  SOLUTIOM 

CALL  RICAT  ( PHI , NPH I , D ,ND, NCONT , KT , NKT , P,NP , DUM ,KDUM) 

CALL  TRANP  ( KT ,NKT , K , NK) 

CALL  PRNT  (P,NP,4HP    ,1) 

CALL  PRNT  (K,NK,4H<    ,1) 
C  CHECK  NOPT  TO  SEE  IF  DONE 

IF  (N0PT.EQ.2)  GO  TO  12 
C 
C 


A  3 


c  ********$**#*#*«****$**$******* 

c  *  * 

C  *   STATE  COVARIANCE  SOLUTION   * 

C  *  * 

£  ******************************* 

C 

C 

READ  (5,14)  (NCGNT(  I )  ,1=1  ,3) 
C  USE  MATRIX  F  FOR    (A-BL),  PHI  FOR  ITS  TRANSPOSE 

CALL  MULT  ( B , NR , L , ML » PH I , NPHI ) 

CALL  SUBT  (  A,  N  A,  PHI  ,  NPHI  ,F,NF) 

CALL  TRANP  ( F , NF , PH I , NPHI ) 
C  SETUP  R  AS  A  ZERO  MATRIX 

CALL  SCALE  ( R , NR , R, NR, O.DOO ) 
C  SET  M=0 

CALL  SCALE  < M , NA , M, NM , O.OOO) 
C  FIND  EXPONENT 

CALL  AUG  ( PHI , NPHI , B , NB, R ,NR, KT , N<T , G,NG, D, NO, Z ,NZ, 0 ) 
C  FIND  EXP 

CALL  ETPHI  (Z,NZ,1.D00,PHIT,NPHIT,DUM,KDJM) 
C  FIND  RICCATI  SOLUTION 

CALL  RICAT  ( PH IT , NPH IT , D, ND,NCONT, Z  1 , NZ1 , M, NM,DUM,KDUM ) 
C  FIND  P*M,  PRINT 

CALL    ADD    (M,NM,P,NP,MP,NMP) 

CALL     PRNT     <MP,NMP,<HP  +  M    ,1) 
C    FIND    CVU,     PRINT 

CALL  TRANP  ( L , NL , LT , NLT ) 

CALL  MULT  (  M,  .MM,  L  T  ,  NLT  ,  MLT  ,  NMLT  ) 

CALL  MULT  ( L , NL , MLT , NMLT, CVU, NC VU) 

CALL    PRNT    (CVU,NCVJ,4HCVJ     tl) 
C    FIND    CVHX,    PRINT 

CALL    MULT    ( h , NH, MP , MMP , PH  I , NPHI  ) 

CALL    MULT     ( PH I ,NPH I , HT ,NHT , PH IT ,NPH IT ) 

CALL     PRNT     ( PHIT,NPHIT,4HCVHX,1 ) 
C  FIND    J 

CALL    TRANP     ( C , NC , CT , NCT ) 

CALL    MULT    ( CT ,NC T, Q, NQ , PH I , NPHI ) 

CALL    MULT    ( PH  I  ,MPHI  , C , NC , PH IT , NPHI T) 

CALL    MULT    ( PH I T , NPHI T , P ,NP , Z , NZ ) 

CALL    MULT     ( S , NS, P, MP , PHI , NPH I ) 

CALL  TPANP  (H,NH,HT,NHT) 

CALL  MULT  ( PHI , NPHI , HT, NHT , PHI T , NPHI T ) 

CALL  MULT  ( PH I T , NPHI T , GI , NGI , PH I , NPHI ) 

CALL  MULT  ( PH I ,MPHI , H , NH, PH I T , NPHI T ) 

CALL  MULT  (PHIT,MPrtI T, P,NP, PHI , NPHI ) 

CALL  TRCE  (Z,NZ,T*A) 

CALL  TRCc  ( PHI  ,MPHI,TRB) 

J  =  TRA+TRB 

CALL  LNCNT  (2) 

WRITE  (6,15)  J 

IF  (N0PT.LT.4)  GO  TO  12 
C 

C 

C  ************************************ 

C  *  * 

C  *   OBTAIN  ROOTS  OF  DET( S I-( A-BL) )   * 

C  *  * 

£  ***** ******* ****** ****************** 

C 

c 

C      SINCE  MATRIX  F  ALREADY  EQUALS  (A-BL)  PROCEED  DIRECTLY 

CALL  DIMCH  (F,NA,DC) 

CALL  CHREQA  ( DC ,NA { 1 ) , CH ) 

CALL  PROOT  (NA(l) ,CH, U,V,1) 

WRITE  (6,16) 

WRITE  (6,17) 

WRITE  (6,18) 

II  =  NA(1) 
C 

DO  3  1=1,11 

WRITE  (6,19)  U(  l),V/(  I  ) 
3  CONTINUE 


A  k 


c 

£  *********************  **  ************* 

r  *  * 

C  *       OBTAIN    ROGTS    OF    DE T( S  I-( A-KH)  )       * 

C  *  * 

Q  ************************************ 

c 

C 

C  FIND    K*H,    PUT     IN    G 

CALL    MULT     ( K, NK, H, NH , G ,N3 ) 
C  FIND    A-KH,     PUT     IN    C 

CALL     SUBT     ( A,NA,G,NG,C,NC) 
C  FINO    ROOTS    AND    P^IMT 

CALL    DIMCH    (C,NC,0C) 

CALL  CHREOA  ( DC, NC ( 1 ) , CH) 

CALL  PROOT  (NC( 1)  ,CH,U,V,1) 

WRITE  (6,16) 

WRITE  (6,23) 

WRITE  (6,13) 

J  J  =  NC(1) 
C 

DO  4  1=1, JJ 

WRITE  (6,19)  U(I),V(I) 
4  CONTINUE 
C 

c 

c  ************************************** 

c  *  * 

C  *   OBTAIN  U(S)/Z(S)  TRANSFER  MATRIX   * 

C  *  * 

Q  ************************************** 

c 
c 

C      FIND  A-BL-KH,  PUT  IN  C 

CALL  SUBT  (F,NA,G,MG,C,NC) 

CALL  DIMCH  (CtNCtDC) 

CALL  CHREQ  (DC ,NC ( 1 ) , CH , 1  ) 

CALL  PROOT  (NC<1)  iC-ltUfVtll 

KK  =  NC(1) 

K2  =  KKH 

WRITE  (6,16) 

WRITE  (6,21 ) 

WRITE  (6,18) 


C 
C 


C 
C 


C 

c 


DO  5  1=1, KK 

5  WRITE  (6,19)  U(I) ,V( I ) 

WRITE  (6,22) 

WRITE  (6,23)  (CH( I) , 1=1, K2) 

LL  =  NL(1) 

DO  11  K1=1,LL 

CAlL  DIMCH  (L,NL,DC) 

MM  =  NL(2) 

DO  6  1=1, MM 

6  W3(I)  =  DC(  I,K1) 

WRITE  (6,16) 
WRITE  (6,24) 
NN  =  NK(2) 

DO    10    Ji=l,NN 

CALL    DIMCH    (K,NK,DC) 

IJ    =    NK(1> 

DO    7    11  =  1, IJ 

7  W4(I I )    =    DC(I  I,J1) 

CALL    MPY    (W4,W3,NK(l) ,CH, NO) 

WRITE    (6,25)     K1,J1 

N2    =    NO  +  1 

CALL    PROOT     (NO,CH,U,V,  I) 


A   5 


c 
c 


10 

11 

12 


13 
14 
15 
16 
17 
18 
19 
20 
21 
22 

23 

24 
25 
26 


WRITE  (6,26) 
WRITE  (6,18) 

DO  8  18=1, NO 

WRITE  (6,191  U(  IS)t\H!8) 

CONTINUE 


DO  9  1=1, N2 
CH(I)  =  -CH( I) 

WRITE  (6,22) 

WRITE  (6,23)  <CH(I) ,I=l,N2) 

CONTINUE 

CONTINUE 
CONTINUE 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
1W-ERS  OF 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
END 
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RDTITL 

SUBROUTINE    POTITl 
COMMON    /LINES/ NLPtLINtTITLE(23) 
READ       (5,100)        (TITLE(I)  ,1  =  1,18) 
100    FORMAT       (18A4) 
CALL    LNCNT(IOO) 
RETURN 
END 
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PRNT 

SUBROUTINE       PRNT ( AR , NAR  ,NAM , I P) 
C         SUBR    PRNT         PRINTS    DOUBLE    PRECISION    MATRIX 

COMMON       /F0RM/NEPR,FMT1(6) , FMT2  16) 

COMMON/LINES/NLP,LIN,TITLE(23) 

COMMON       /MAX/MAXRC 
C-    NOTE    NLP    NO.     LINES/PAGE    VARIES    WITH    THE     INSTALLATION. 

DATA    KZ,KW,KB    /1H0,IH1,1H    / 

REAL*8    AR 

DIMENSION    AR(1 ) ,NAR(2) 

NAME  =  NAM 
C-IF  IP  =1, HEADLINE  SAME  PAGE,  IF  IP  =2,   HEADLINE,  NEW  PAGE 
C      IP=3t   NO  HEADLINE,  SAME  PAGE,     IP=4,  NO  HEADLINE,  NEW  PAGE 

I  I  =  IP 

NR=NAR(1) 

NC=NAR<2) 

NLST  =  NR  *  NC 

IF(NLST  .GT.  MAXRC  .OR.  NLST  .LT.  1 . OR.NR. LT . 1 )   GO  TO  16 

IF1NAME   .EG.  0)  NAME  =  KB 
C-  SKIP  HEADLINE  IF  REQUESTED. 

GO  TO    (11,10,132,12) ,        II 

10  CALL  LNCNT( 100) 

11  CALL  LNCNT( 2) 

3  WRITE (6, 177)  KZ,  NAME , NR, NC 
177    FGRMAT(A1,5X, A4,8H   MATRI X , 5X , 13 , 5H  ROWS ,5X , I  3 , 3H  COLUMNS ) 
GO  TO  13 

12  CALL  LNCNT( 100) 
GO  TO  13 

132  CALL  LNCNT(2) 
WRITE  (6,891) 
891  FORMAT  (1H0) 
C-  BELOW  COMPUTE  NR  OF  LINES/  ROW  —DECIDE  IF  1  EXTRA  BLANK  LINE 

13  J=(NC-1 )/NEPR*l 

C-  WHY  ALWAYS  ADD  1  LINE-  BECAUSE  IF  MULTIPLE,  USE  1  BLK  LINE  EXTRA. 

NLPW  =  J 

JST  =  1 
C-  COMPUTE  LAST  ROW  P3SITI0N  -1  BELOW 

NLST  =  NLST  -NR 

MU=NC 

IF   (NC.GT.NEPR)    *1N  =  NEPR 

KLST=NR*(MN-1) 
91     CONTINUE 

DO  912  J  =  JST,  NR 

CALL  LNCNT(NLPW) 

KLST  =  KLST  +1 

WRITE  (6,FMT1)     (AR(N),  N  =  J,  KLSTt  NR ) 

IF   (NC.LE.NEPR)    GO  TO  912 

NLST  =  NLST  +1 

KNR=KLST+NR 

WRITE  (6,FMT2)(AR(N) , N=KNR , NLST , NR J 
912    CONTINUE 

RETURN 
16  CALL  LNCNT( 1) 

WRITE   (6,916)    NAM, NAR 
916  FORMAT   <■  EPROR  IN  PRNT   MATRIX  ',A4,»  HAS  NA=*,2I6) 

CALL  ASPERR 

RETURN 

END 
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LNCNT 

SUBROUTINE    LNCNT       (Ni 
COMMON/LI  NFS/ NIP, LIN, TITLE  (23) 
L  IN=LIN+N 

IF       <  LIN.LF.MLP)  GO    TO    20 

WRITE       (6,1010)       ( TITLE( I )  ,1  =  1,23) 
1010    FORMAT       (  1H1,23A4/) 
LIN=2+N 

IF       (N.GT.NLP)       LIN  =  2 
20    RETURN 
END 
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300 

999 

50 

1000 


SUBROUTINE  ADD 
DIMENSION  A(l) ,»(1) 
COMMUN   /MAX/MAXRC 
DOUBLE  PRECISION  A, 
IF( (NA{ 1) .NE.N3I 1 ) ) 
NC< 1)=NA(1) 
NC(2) =NA(2) 
L=NA(  1)*NA( 2) 
IF   (NA(l).LT. 
DO  300  I- 1 , L 
C(I)  =  A( I)+B(  I  ) 
GO  TO  1000 
CALL  LNCNT  (1) 
WRITE(6,50)  NA 


ADD 

< A,NA.B,NB,C.NC) 

,C( 1),NA(2) ,NB(2) ,NC<2) 


R  C 

.CR.(NA(2) .NE.N8(2)) )  GD  TO  999 


l.OR.L.LT.l.OR.L.GT.MAXRC) 


10  TO  999 


■  NB 


FORMAT   (•  DIMENSION  ERROR  IN  ADD 

CALL  ASPERR 

RETURN 

END 


NA=« ,2l6t5X,f NB=»,2I6) 
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Main: 

SUBROUTINE  MULT  (  A  ,  NA  ,  B  ,  NB  ,  C  ,  NO 

DIMENSION  A(l),6(l),C(I),NA(2),NB(2),NC(2) 

DOUBLE  PRECISION  A,6,C 

COMMON   /MAX/MAXRC 

NC<1)  =  NA(  1  ) 

MC(2)  =  NB(2) 

IF(NA(2).NE  .NB< 1) )  30  TO  999 

NAR  =  NA(l ) 

NAC  =  NA(2) 

NBC  =  Nb(2) 

NAA*NAR*NAC 

NBB=NAR*NBC 

IF   (NAR.LT.1.GR.MAA.LT.1.0R.NAA.3T.MAXRC.0R.NBB.LT.1 .OR. 
1     NBB.GT.MAXRC)   GO  TO  999 

IR  =  0 

IK=-NAC 

DO  300  K=1,NBC 

IK  =  IK  +  NAC 

DO  300  J=1,NAR 

I R=I R>1 

IB=IK 

JI=J-NAR 

C(IR)=0. 

DO  300  1=1, NAC 

J I  =  J I  +  NAR 

I  3=1 B-H 

C( IR) =C( IR) +A( JI )*B( IB) 
300  CONTINUE 

GO  TO  1000 
999  CALL  LNCNT  (1  ) 

WRITE (6, 500)  <NA( I) , 1=1,2)  ,<NB(  I) ,1  =  1,2) 
500  FORMAT   (•  DIMENSION  ERROR  IN  MULT     NA=« , 21 6, 5X , • NB=« t 21 6 ) 

CALL  ASPERR 
1000  RETURN 

END 
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SCALE 

SUdROUTIME  SCALE  (A,  NA.  B,  N3,  S) 

DIMENSION  A(l) ,3(1) ,NA(2) ,NB(2) 

COMMON   /MAX/MAXRC 

DOUBLE  PREC  ISION  A*  B,  S 

N8(l )  =  NA(  1) 

NB(2)  =NA(2) 

L  =  NA( 1)*MA(2) 

IE   (NAU ) .LT.l.OR.L.LT.l.OR.L.GT.MAXRC)   30  TO  999 

DO  300  1=1, L 
300  P (I)=A( I)*S 
1000  RETURN 
999  CALL  LNCNT( 1J 

WRITE   (6,50)  NA 
50  FORMAT   (*  DIMENSION  ERROR  IN  SCALE    NA=«,2I6) 

CALL  ASPERR 

RETURN 

END 
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TRANP 

SUBROUTINE  TRANP ( A . NA T 6.NB) 

DIMENSION  A<1) ,B(1) ,NA(2) ,NB<2) 

DOUBLE  PRECISION  A,B 

COMMON   /MAX/MAXRC 

N8(1)=NA(2) 

NB(2) =NA( I) 

NR=NA( 1) 

NC=NA(?) 

L=NR*NC 

IF   (NR    .LT.1.0R.L.LT.L.0R.L.GT.MAXRC)   30  TO  999 

IR=0 

DO  300  1=1, NR 

IJ=I-NR 

DO  300  J=1,NC 

I J=IJ+NR 

IR=IR+1 
300  B(IR)=A(IJ) 

RETURN 
999  CALL  LNCNT(  1) 

WRITE   (6,50)   NA 
50  FORMAT   <•  DIMENSION  ERROR  IN  TRANP    NA=»,2I6) 

CALL  ASPERR 

RETURN 

END 
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SUBROUTINE     INV ( A , NA . DET  .  L ) 
DIMENSION    AU),    L(ll,NA(2) 


DOUBLE    PRECISION    A,     DET,    BIGA    ,       HOLD 
COMMON       /MAX/MAXRC 
IF       (NAUI.NF.NA12))        GO    TO    600 
C  SEARCH    FOR    LARGEST    ELEMENT 

D£T=  1. 
N=NA( 1) 
NSQ=N*N 

IF       (N.LT.l .OR.NSQ.GT.MAXRC)  GO    TO    600 

NK    =    -    N 
DO    80    K=    1,    N 
NK    =    NK    ♦    N 
LIK)     =    K 
NPK=N+K 
L(NPK)=K 
KK    =    NK    +    K 
BIGA    =    A(KK) 
DO    20    J=    K,     N 
IZ    =    N*(J    -     1) 
DO         20    1=    K,    N 

u  =  iz  +  r 

10    IFOABS(BIGA)     -    DA3S(A(IJ1))     15,    20,     20 
15    BIGA    =    A(  IJ  ) 

L(K)     =    I 

NPK=N+K 

L(NPK)=J 
20  CONTINUE 
C      INTERCHANGE  ROWS 

J  =  L(K) 

IF(J  -  K)  35,  35,  25 
25  KI  =  K  -  N 

DO  30  I  =  1,  N 

KI  =  KI  ♦  N 

HOLD  =  -A(KI) 

JI  =  KI  -  K  +  J 

MKI  )  =  A(J1I 
30  A( JI )  =  HOLD 
C      INTERCHANGE  COLUMNS 
35  NPK=N+K 

I=L(NPK) 

IF  (I  -  K)  45,  45,  33 
3  8  J  P  =  N*  (  I  -  1  ) 

DO  40  J=  1,  N 

JK  =  NK  +  J 

JI  =  JP  +  J 

HOLD  =  -A( JK) 

A( JK)  =  A( J  I) 
40  A(  JI  )  =  HOLD 
C      DIVIDE  COLUMN  BY  MINUS  PIVOT(VALUE  OF  PIVOT  ELEMENTS  IS  CONTAINS 
C      IN  BIGA) 

45  IF(BIGA)  48,  4b,  48 

46  DET  =  0. 

CALL  LNCNT  (1) 

KKK=KK-1 

WRITE   (6,1046)   KKK 
1046  FORMAT   (•  INV  ERROR   DETERMINANT  OF  A=0    RANK  OF  A=*,I4) 

CALL  ASPERR 

RETURN 
48  DO  55  1=  1,  N 

IF  ( I  -  K)  50,  55,  50 
50  IK  =  NK  *  I 

AUK)  =-AUK)/(BIGA) 
55  CONTINUE 
C      REDUCE  MATRIX 

DO  65  1=  1,  N 

IK  =  NK  +  I 

HOLD  =  AUK) 

I  J  =  I  -  N 

DO  65  J=  1,  N 

IJ  =  IJ  +  N 

I F ( I   -  K  )  60,  65,  60 
60  IF(J-  K)  62,  65,  62 
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INV   (cont'd) 

62    KJ    =     IJ    -I     ♦    K 

A(  IJ)     =    HOLD*    A(KJ)     ♦    A(  I  J) 
65    CONTINUE 

DIVIDE    ROW    BY    PIVOT 

KJ    =    K    -    N 

DO    75    J=    1»     N 

KJ    =    KJ    *    N 

IF(J  -  K)  70,  75,  70 
70  A(KJ)  =  A(KJ)/BIGA 
75  CONTINUE 

PRODUCT  OF  PIVOTS 

DET=DET*BIGA 

REPLACE  PIVOT  BY  RECIPROCAL 

A(KK)  =  l./BIGA 
80  CONTINUE 

FINAL  ROW  AND  COLUMN  INTERCHANGE 

K  =  N 
100  K  =  K  -  1 

IF(K)  150,  150,  105 
105  I  =  L(K) 

IF  (I  -  K  )  120,  120,  108 
108  JQ  =  N*(  K  -  1) 

JR  =  N*( I-  1) 

DO  110  J=  1,  N 

JK  =  JQ  ♦  J 

HOLD  =  AUK) 

JI  =  JR  ♦  J 

A( JK)  =  -  A(JI  ) 
110  A( JI )  =  HOLD 
120  NPK=N+K 

J=L(NPK) 

IF(J  -  K)  100,  100,  125 
125  KI  =  K  -  N 

DO  130  1=  1  ,  N 

KI  =  KI  +  N 

HOLD  =  A(KI  ) 

JI  =  KI  -  K  +  J 

A(KI )  =  -  A(JI  ) 
130  A(JI)  =  HOLD 

GO  TO  100 
150  RETURN 
600  CALL  LNCNT  (1) 

WRITE   (6,1630)   NA 
1600  FORMAT   <•  INV  REQUIRES  SQUARE  MATRIX    NA=«,2I4) 

CALL  ASPERR 

RETURN 

END 
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NORM 

SUBROUTINE  NOPM( A. NA , ANORM ) 

DIMrNSION  NA(2)  ,A(  1) 

DOUBLE  PRECISION  A , ANORM, SUM, ROWMAX ,COLMAX 

CUMMON   /MAX/MAXRC 

NAR  =  NA<1) 

NAC  =  NA(2) 

L=NAR*NAC 

IF   (NAR   . LT.l.OR.L.LT.l .OR.L.GT.MAXRC)   GO  TO  999 

COLMAX  =  0. 

ROWMAX  =  0. 

K  =  0 

DO  300  I  =  l.NAC 

SUM  =  0. 

DO  301  J  =  1,NAR 

K  =  K  *  1 

301  SUM  =  SUM  +  PA3S(A(K)  ) 

IF  (COLMAX.LT.SUM)  COLMAX  =  SUM 
300  CONTINUE 

DO  302  I  =  l.NAR 

SUM  =  0. 

K  =  I  -  MAR 

DO  303  J  =  l.NAC 

K  =  K  *  NAR 
303  SUM  =  SUM  +  DA3S(A(K)) 

IF    (ROWMAX. LT. SUM)     ROWMAX=SUM 

302  CONTINUE 

ANCRM    =    DMINKCCLMAX,  ROWMAX) 
RETURN 
999    CALL    LNCNT     (1) 

WRITE       (6,50)       NA 
50    FORMAT       («     DIMENSION    ERROR    IN    NORM  NA=*,2I6) 

CALL    ASPERR 
RETURN 
END 
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UNITY 

SUBROUTINE    UN  I TY ( A , NA ) 

DIMENSION    A( 1) ,NA( 2) 

DOUBLE  PRECISION  I 

IF(NAU)  .NF.NA(2)  )  GO  TO  999 

CALL  SCALE(  A,N A, A, NA , O.DO ) 

J  =  -  NA( 1) 

NAX  =  NA(1) 

DO  300  1=1, NAX 

J=NAX  +  J+1 
300  A( J)=l. 

GO  TO  1000 
999  CALL  LNCNT  (1 ) 

WRITE16,  50)(NA( I), 1=1,2) 
50  FORMAT   (•  DIMENSION  ERROR  IN  UNITY    NA=»,2I6) 

CALL  ASPFRR 
1000  RETURN 

END 
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EQUATE 

SUBROUTINE     EOUAT E< A, NA . 3 . NB ) 

DIMENSION    All)  ,B(1) ,NA(2J  fNB(2) 

DOUBLE  PREC  ISIGN  A,  B 

COMMON   /MAX/MAXRC 

NB(1)  =  NA(  1) 

NB(2)  =NA(2) 

L=NA( 1)*NA(2) 

IF   (NA( 1 ). LT. l.OR.L.LT.l.OR.L.GT.MAXRC)   SO  TO  999 

DO  300  1=1, L 
300  B(I)=A( I) 
1000  RETURN 
999  CALL  LNCNT  (1) 

WRITE   (6,50)   NA 
50  FORMAT   (•  DIMENSION  ERROR  IN  EQUATE   NA=»,2I6) 

CALL  ASPERR 

RETURN 

END 
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ETPHI 

SUBROUTINE     ETPHI ( A , NA , TT, B ,N8 t DUMMY , KDUM) 
01  MENS  ION    A( 1 )  ,B(1)  , DUMMY*  I ) , NA ( 2 ) f MB ( 2 )  , ND ( 2 ) 
DOUBLE    PRECISION    A,     T,    TT,     ANAA,    TMAX,    8,    DUMMY    ,S,     SC 
COMMON       /MAX/MAXRC 
NR=NA(1) 
NCC=NA(2) 
NBU)=NR 
NB(2) =  NCC 
LD=NP*NCC 

IF     (NR.NE.NCC.OR.NJR.LT. 1  .0R.LD.3T.MAXRC)     GO    TO    998 

NDD=2*NA(1)  *NA(i ) 
IFUDUM       .LT.NDD)     GO    TO    998 
NDD=    NA(1)*NA( 1)+1 
T  =  TT 

CALL    NORM(A,NA,ANAA> 
TMAX=     10.01/ANAA 
K  =  0 
iOi     IF    (TMAX-T     )     103,104,104 

103  K  =  K+1 
T=TT/2**K 

IF    (K-1000)     101,102,102 

104  SC=T 

CALL     SCALEC A,NA,A,NA,T) 

CALL    UNITY(BfNB) 

11=2 

N    =    35 

CALL     ADD(A,NA,B,NB,DUMMY(1 ),N0) 

CALL    EQUATE(A,NA,DUMMY(NDD) ,ND) 

106  CALL    MULT(A,NA,DUMMY(NDD) ,N0,3,NB) 
SM.DO/I  I 

CALL  SCALE(B,NB,DJMMY(NDD) ,ND,S) 

CALL  ADD(0'JMMY(\IDD)  ,ND,DJMMY(  I)  ,ND,B,NB) 

CALL  EQUATE(S,NB,DUMMY(1),ND) 

N  =  N-1 

IF  (N)  107,107,105 

105  11=11+1 
GO    TO    106 

107  IF    (K)     109, 108,212 

109  CALL     LNCNT    (1  ) 
WRITE    (6,110) 

110  FORMAT    (•     ERROR    IN    ETPHI       K    IS    NEGATIVE') 
112    IF    (K-l)     213,212,212 

213    K=l 

212    DO    11  1    J=1,K 

T  =  2*T 

CALL    EQUATE(3,    NB ,    DUMMY(l),    ND ) 

CALL     EOUATE(OUN<MY(  1  )  ,    ND,     DUMMY(NDD),    ND ) 

111  CALL    MULT(DUMMYIIMDD)  ,ND,DUMMY(  1  )  ,ND,B,NB) 
TT  =  T 

108  CONTINUE 
S=1.D0/SC 

CALL     SCALE( A,NA,A,NA,S) 

RETURN 
102    CALL     LNCNT     (1) 

WRITE    (6,119) 
119    FORMAT  ('     ERROR    IN    ETPHI  K=1000») 

CALL    ASPERR 

RETURN 
998    CALL     LNCNT    (I ) 

WRITE       (6,50)       NA,KDUM,NDD 
50    FORMAT       (•     DIMENSION    ERROR    IN    ETPHI  NA=',2I6,  • KDUM= • , I  5 , 5X, 

1  •KDUM(MIN}=» , 15) 

CALL    ASPERR 

RETURN 

END 
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AUG 

SUBROUTINE    AUG< F ,NF , G, NG , RI , NRI ,H,     NH.Q,NQ,     C,NC,Z,NZ,     II) 
DIMENSION    F(  1  )  ,G( 1) ,^1  (  1)  ,H( 1)  ,Q(  1 )  ,Z(  1)  ,C(  1  ) 

DIMENSION    NNP1(2),NNP2(2) , NNP3(  2 ) , NNP4( 2 ) , NF ( ? )  ,NG ( 2 )  , NRI ( 2  ) , 
1NH<2) ,NZ(2) ,NC<2) ,NN<2) ,NQ(2) 
DOUBLE    PRECISION    t ,     G,    RI,H,Q,C,Z 
IF((NF(1) ,NE.NF<2) ).OR. (NRI (1 ).NE.NRI<2) ).0R.( 
INGU  )  .NE.NQ(2)  )  )     GO    TO    995 
NNZ=2*NF(1) 

IF(     (NG(1) . NE.NF< 1) ) . OR. ( NG ( 2 ) . NE . NR I ( 1 ) ) )30    TO    995 
IF(  II  .EQ.U     GO    TO    206 

IF((NH< 1 ) .NE.NQil )) .OR.(NH(2) .NE.NF(2)) )    G3    TO    995 
206    TWO    =    2 
N    =    NF( 1) 
NSQ    =    N*N 
NZ(1)=NNZ 
NZ(2)=NNZ 
NP1  =  1 

NP2    =    NP1     f    NSQ 
NP3    =    NP2+NSQ 
NP4    =    NP3    +    NSQ 
CALL     TRANP(G,NG,Z(NP4) ,NNP4) 
CALL    MULT (R  I, NRI  , Z(NP4) ,NNP4,C,NC) 
CALL    MULT(G,NG,C,NC,Z(NP4) ,    NNP4) 
IF( I  I     .EQ.     1)    GO    TO    204 
CALL     TRANP<  H,NH,Z(,\|P3)  ,    NNP3) 
CALL    MULT( Q,NQ,H,NH,Zl 1) ,     NNP1 ) 

CALL    MULT( Z(NP3) ,NNP3,Z(       1 ) ,NNP1 , Z < NP2) ,NNP2) 
GO    TO    205 

204  CALL    EQUATEIQ,    NQ,    Z(NP2),     NQ) 

205  NPAIR=M0D(N,2) 
IF(NPAIR.EQ.O)     NPAIR=TWO 
NN( 1 )     =    N 

NN(2)     =    1 

GO    TO    (201,202),NPAIR 

201  DO    300    1=1, N, 2 
NP2    =    N*(N+I-1 )+l 

N  TH3=  TWO*N*  ( I -1 ) +N+ 1 

300  CALL    EQUATE (Z(NP2) , NN , Z ( NTH3 ) , NN) 
DO    30  2    1=2, N, 2 
NP4=N*(3*N*I-1 )  +  l 
NTH2=TW0*N*(N+I-1)+1 

302  CALL    EQUATE (ZCNP4) ,NN,Z(NTH2) ,NN) 
GO    TO    ( 202, 203) ,NPAIR 

202  DO    301    1=2       ,N,2 
NP2    =    N*(N+I-1)<-1 
NTH3=TW0*N*< 1-1 )*N>1 

301  CALL    EQUATE(Z(NP2) , NN, Z ( NTH3 ) , NN) 
DO    304    1=1, N, 2 
NP4=N*(3*N+I-1)+1 
NTH2=TW0*N*(N+I-1)+1 

304    CALL    EQUATE(Z(NP4) ,NN,Z(NTH2) ,NN) 
GO    TO     (203,201) ,NPAIR 

203  DO    303    1=1 ,N 
I  J    =     I+N 

DO    303    J=1,N 
JJ    =    J+N 
L1=NNZ*( J-l )+I 
L2=NN2*( I  J- II *JJ 
L3=N*( J-l)fI 
Z(L1 )=-F(L3) 

303  Z(L2)=F(L3) 
GO    TO    1000 

995    CALL     LNCNT     (2) 

WRITE       (6,50)       NF,NG,NRI,NH,NQ 
50    FORMAT       (•     DIMENSION    ERROR    IN    AUG' , T 35 , ' NF • , 9X , • NG1 , 9X , «NRI • , 8X, 
1  'MH' ,9X,  •NQV29X,  5(3X,2I6)  ) 

999    CALL     ASPERP 
1000    RETURN 
END 


A  20 


RICAT 

SUBROUTINE    R IC AT ( PH I T NPH I .C ,NC , NCONT ,K, NK , PT ,NPT, W, KOUM ) 
01  MENS  I  ON    NCONT(i) ,NPhI (2) ,NC(2) ,NK (2),NN(2) ,NM(2),NPT(2) 
DIMENSION    PHI  ( 1)  ,C( 1 )  ,K(  1)  ,PT<  1 ) TW( 1  ) 

DOUBLE    PRECISION    PHI,     C,    K,    PT,     SUM,     SUMN,    AL,     W,DET 
TW0=2 

N    =    NPHK  D/TWO 
NSQ=N*N 
NSQ4=4*NSQ 
NP1=1 

NP2=    NSQ+NP1 
NP3=NSQ+NP2 
NP4=    NSQ+NP3 

IF       (KDUM.LT.NSy4     )        GO    TO    600 

IF   (NPH I (2) .N E.N  PHI ( 1 ).0R.NC(2) .NE . N.OR.NPT ( 1 ) .NE. N. OR.NPT ( 2) 
1   NE.N)    GO  TO  600 
NPRINT=NCONT( 1 ) 
N8LOCK=NCCNT(2)/NPRINT 
NZ=NC0NT(3) 
C      REARRANGE  PHI  MATRIX 
NN(1) =N 
NN(2)=1 
DO  300  1  =  1, N 
NTH1  =TWO*N*( I-l)+l 
NTH3=NTH1+N 
NW1  =  IM*( 1-1) +1 
NW2=NW1+N*N 
CALL  EQUATE ( PHI ( NTH  I ) , NN, W( NW1 )  ,NN ) 

300  CALL  EQUATE(PHI(NTH3) ,NN,W(NW2)  ,NN) 
NM(1)=TW0*N*N 

NM(2)=1 

CALL  EQUATE(W( I) ,MM,PHI( 1) ,NM) 

DO  301  1=1, M 

NTH2=TW0*N*(N+I-l)f 1 

NTH4=NTH2+N 

NW3    =    N*(TW0*N+I-1 )+l 

NW4=    NW3+N*N 

CALL    EQUATE (PHI (NTH 2 ) ,NN,W(NW3) ,NN) 

301  CALL     EQUATE(  PHMNTH4  )  ,NN,W(NW4)  ,NN) 
NWW=Twu*N*N+i 

CALL     EQUATE ( W( NWW ) , NM, PHI (NWW) ,NM) 
C  MAIN    COMPUTATION 

C  CALL    UNITYt PT,NPT) 

DO    306    1=    1  ,N 
306    K(I)     =    0. 

NTOT=0 

DO    403    I=1,NBL0CK 

DO    104    J=1,NPRINT 

CALL  MULT(PHI(NP3) ,  MPT,  PT,  NPT,  W(l),  NPT) 

CALL  ADD  (PHI(l),  NPT,  W(l),  NPT,  W(l),  NPT) 

CALL  IiMV(W(l},  MPT,  DET,  W(NP2)) 

CALL    MULT(PHI (NP4) ,     NPT,    PT,    NPT,     W(NP2),     MPT) 

CALL    ADD(PHI(NP2)  ,    NIPT,    W(NP2),    NPT,    W(NP2),     NPT) 

CALL    MULT(W(NP2),    NPT,    W(l),    NPT,    PT,    NPT) 

SUMN=0. 

SUM=0. 

DO    202    IL  =  1  ,N 

ILP=IL+NP3 

NIL=(  IL-1)*N>IL 

SUM=SUM+DABS(PT(NIL> ) 
202    SUMN=SUMN*DABS(PT(NIL)  -W(ILP)) 

AL=SUMN/SUM 

DO    201     IL  =  1  ,N 

NIL  =  (  IL-1  )*N+IL 

ILP=IL+NP3 
201  W(ILP)    =PT(NIL) 
204  DO  104  M=2,N 

N1=M-1 

DO  104  L=1,N1 

L1=N*(L-1)  *M 

L2=N*(M-1)+L 

PT(L1)=(PT( LI)  +  PT(L2)  )/2. 

PT(L2)=PT(L1) 

IF(AL-. 00001)  203,203,104 
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RICAT    (cont'd) 

104    CONTINUE 

NTOT=I*NPRINT 

GO    TO    305 
203    NTOT=NTOT+J 

305    CALL       MULT     ( C , NC , PT , NPT ,K,NK) 
103    GO    TO     (404,400,401,402) t     NZ 

400  CALL     LNCNT     ( 1) 
WRITE     (6,     50)     NTOT 

50    FORMAT    (10X,I4,14H  ITERATIONS       ) 

CALL     PRNT        <PT,NPT,' P( T) •  ,1) 
GO    TO    403 

401  CALL    LNCNT     (1) 
WRITE    (6,     50)    NTOT 
CALL     PRNT  (K,NK,f K(T)' ,1) 
GO    TO    403 

402  CALL    LNCNT     (1  ) 
WRITE     (6,     50)     NTOT 
CALL    PRNT  (K.NKf'KJTIM) 
CALL    PRNT       < PT,NPT,'P(T)  •  ,1  ) 
IFCAL-. 00001)     210,210,403 

404  IFCAL-. 00001)     ^05,405,403 

403  CONTINUE 

405  CALL     LNCNT     ( 1 ) 
WRITE(6,50)NT0T 
REARRANGE    PHI    MATRIX 

210    CALL     EQUATE(PHI( 1)  ,MM,W( 1  )  ,NM) 
DO    303    1=1 ,N 
NTH1    =TWO*N*( I-l)+l 
NTH3=NTH1+N 
NW1=N*<  I-D+l 
NW2  =  NWH-N*N 
CALL     EQUATE ( W( NW1 ) , N N ,PHI (NTH1 ) ,NN) 

303  CALL    EQUATE  (  W(NW2)  ,NN,PHI  ( NTH3 )  ,M;M) 
CALL     EQUATE  (  PHHNWW)  ,NM,  W(NWW)  ,NM) 
DO    304    1  =  1, N 
NTH2=TWQ*N*(N*I-1)+1 
NTH4=NTH2+N 
NW3    =    N*(TW0*N4-l-l  )*1 
NW4=    NW3*N*N 
CALL    EQUATE( W(NW3) ,NN,PHI (NTH2) ,NN) 

304  CALL    EQUATE! W(NW4) ,NN , PH I ( NTH4 ) ,NN) 
RETURN 

600    CALL    LNCNT     (2) 

WRITE       (6,1600)       NPHI , NC , NPT , KDUM, NS04 
1600    FORMAT       ('     DIMENSION    ERROR    IM    RICAT •  ,T35 , ' ^PHI •  ,  7X, • NC « ,9X, » NPT» 
1  ,6X, 'KDUM' ,3X, • KDUM (M IN) • /29X , 3 ( 3X , 2 14 ) , 4X , I  4, 5X , 1 4 ) 

CALL    ASPERR 
RETURN 
END 
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ASPERR 

SUBROUTINE  ASPERR 

OATA  I  /10/ 
C      CALL  TRACE 

C     ERRTRA  IS  THE  360/67  TRACE  ROUTINE      TRACE  IS  FOR  TSS 
C      CALL  ERRTRA 

C      THIS  IS  AN  INSTALLATION  DEPENDENT  SUBROUTINE 

C      SUBROUTINE  ERRTRA  IS  A  SUBROUTINE  SUPPLIED  BY  THE  AMES  OPERATING 
C  SYSTEM  TO  PROVIDE  AN  ERROR  WALKBACK 

C      THE  STATEMENT   "CALL  ERRTRA"  SHOULD  BE  EITHER 
C  1)  CHANGED  TO  MATCH  THE  USERS  OPERATING  SYSTEM, 

C        OR  2)  DELETED  ALTOGETHER. 

1=1-1 

IF  (I  .GT.O)   RETURN 

1  =  10 

WRITE  (6,100) 
100  FORMAT  (•  TOO  MANY  ERRORS.    EXIT  CALLED1) 

CALL  EXIT 

RETURN 

END 
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BLOCK  DATA 

BLOCK    DATA 

COMMON      /FORM/NFPR,FMTl (6) f FMT?<6) 
COMMON/LINES/NLP,LINtTITLE<23) 
COMMON       /MAX/MAXRC 
DATA       MAXRC/C4-00/ 
C-    MOTE    NLP    NO.    LINES/PA3E    VARIES    WITH    THE     INSTALLATION. 
DATA    LIN,NLP/1 ,7^/ 

DATA       NEPRiFMTl  / 7 , • ( 1 P7 • , ■ D16 . • , • 7 ) • / 

DATA       FMT2/ «(3X,', •iP7D« , • 16.7' , •) • / 
DATA       TITLE       /19*«  S'OSPAN'C    PR  •  ,  •  OGRA  ■  ,  ■  M 

END 
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READ1 

SUBROUTINE    READ1     ( A, NA, NZ ,NAM) 

COMMON       /MAX/MAXRC 

DIMENSION       A(l        ) , NA ( 2 ) ,NZ < 2 ) 

DOUBLE    PREC ISION    A 

IF       (NZ(  D.EQ.O)       30    TO    410 

NR=NZ(1) 

NC=NZ<2) 

NLST=NR*NC 

IF(NLST    .GT.    MAXRC     .OR.    NLST    .LT.    1 .OR.NR.LT . 1 )       GO    TO    16 

DO    400    I    =     1,     NR 
400    READ    (5tl01)     (A(       J),     J    =    I,NLST,NR) 

NA(1 )=NR 

NA(2)=NC 
410  CALL   PRNT  (A,NA,NAM,1) 
101  FORMAT  (7E10.5) 

RETURN 
16  CALL  LNCNTd  ) 

WRITE   (6t916)   NAM,NR,NC 
916  FORMAT   (•  ERROR  IN  READ1    MATRIX  » ,A4,»  HAS  NA=»,2I6) 

CALL  ASPERR 

RETURN 

END 
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READ 


SUBROUTINE  READ  (I,  A,  NA. 
DIMENSION  A(l ) ,  8(1),  C ( i ) 
DIMENSION  NA(2),  NB ( 2  )  ,  NCX { 2 ) t 
DOUBLE  PRECISION  A,  8,  C,  D,  E 


Bt  NB, 
0(1), 
ND( 


100 

999 


READ< 

CALL 

IF  (  I 

READ( 

CALL 

IF  (  I 

READ< 

CALL 

I  F(  I 

READ( 

CALL 

IF(I 

READ! 

CALL 


NZ 
TO  999 


LAB) 


,NZ,   LAB) 
TO  999 


,100)  LAB, 
READKA,  NA 

EO.  I)  GO 
5,100)  LAB, 
READKB,  NB 

EQ.  2)  GO 

,100)  LAB, 
READKC,  NCX,NZ,  LAB) 
.EQ.  ?)  GO  TO  999 
5,100)  LAB, 
READKD,  ND,NZ,   LAB) 
.EO.  4)  GO  TO  999 
5,100)  LAB, 
PEADKE,  NE,MZ,   LA8) 
F0RMAT(A4,4X,2I4) 
RETURN 
END 


NZ(  1 


NZ(  1 


NZ(  1 


NZ(  1 


NZ(  1 


C,  NCX,  D,  NO,  E,  NE) 

Ed) 

2)  ,  NE(2) ,NZ<2) 

NZ<2) 


NZ12) 


NZ(2) 


NZ(2) 


NZ(2) 
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SUBT 
SUBROUTINE  SUBT(A,.MATBfN9tCfNC» 


(1) 


,t( 1),NA(2),NB(2),NC(2) 
B,C 


OR.(NA(2) .NE.NB(2)))  GO  TO  999 


300 

999 

50 

1000 


DIMENSION  A(1),B 

DOUBLE  PRECISION  A, 

COMMON   /MAX/MAXRC 

IFUNAU ) ,NE.NB( 1 )) 

NC(1)=NA(1) 

NC(2)=NA(2) 

L=NA( 1)*NA(2) 

IF   (NA(l).LT.l.OR.L.LT.l.OR.L.GT.MAXRC)   30  TO  999 

DO  300  1=1, L 

C  (I)=A(I)-B( I) 

GO  TO  1000 

CALL  LNCNT  (1 ) 

WRITE(6,50)  NA 


NB 
FORMAT   («  DIMENSION  ERROR  IN 
CALL  ASPERR 
RETURN 
END 


SUBT 


NA=« ,2I6,5X,«NB=» ,216) 
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(A,NA,TR) 
A( 1 ) ,TR 


GO  TO  600 


10 

600 

1600 


TRCE 

SUBROUTINE  TRCE 

DOUBLE  PREC ISION 

DIMENSION  NA(2) 

COMMON   /MAX/MAXRC 

IF  (NA<1> .NE.NAl  2)  ) 

TR=0.00 

N=NA( 1) 

NN=N*N 

IF   (N.LT.l .OR.NN.GT.MAXRC) 

DO  10  1=1, N 

M=I+N*( 1-1) 

TR=TR+A(M) 

RETURN 

CALL  LNCNT1 1) 

WRITE  (6,1600)  NA 

FORMAT  (»  TRACE  REQUIRES  SQUARE 

CALL  ASPERR 

RETURN 

END 


GO  TO  600 


MATRIX 


NA=,,2I6) 


A  28 


PROOT 

SUBPOUTINEPROOT(N,A,U,V,IR) 
C       THIS    SUBROUTINE    USES    A    MODIFIED    BARSTOW    METHOD    TO    FIND    THE 
C  ROOTS    OF    A     POLYNOMIAL. 

DOUBLE    PRECISION    A ( 20  )  ,U( 20 ) ,V ( 20 ) ? H( 21 ) , B < 2 1) , C( 21 ) , P ,Q ,R , F , E 
1  CBAR,D,QP,PP,ZZ 

DO    91     1  =  1  ,  iN 

U(  I)=0. 

V(  I)=0. 
91    CONTINUE 

IREV=IR 

NC=N+1 

D01I=i,NC 
i    H(I)=A( I) 

ZZ  =  0. 

DO    90    1=2, NC 

ZZ=DAeS<H( I )) fZZ 
90    CONTINUE 

IFUZ.LT.  l.D-10)    30    TO    100 

P  =  0. 

Q=0. 

R=0. 

3  IF(H(1))4,2,4 
2    NC=NC-1 

V(NC)=0. 
U(NC)=0. 
0010021=1, MC 
1002    H( I)=H(H-1) 
G0T03 

4  IF(NC-1J5,100,5 

5  IF(NC-2)7,6,7 

6  R=-H(l)/H(2) 
GOT050 

7  IF(NC-3)9,8,9 

8  P=H(2)/H(3) 
Q=H(1)/H(3) 
G0T070 

9  IF(DABS(H(NC-1)/H(NC) ) -DABS< H< 2  )/H ( 1)))10,19,19 

10  IREV=-IREV 
M=NC/2 
D011I=1,M 
NL=NC+1-I 
F=H(NL) 
H(NL)=H(  I) 

11  H(I)=F 

IF(Q) 13,12,  13 

12  P=0. 
G0T015 

13  P=P/C 
Q=l./G 

15  IF(R)  16,19, 16 

16  R=l./R 

19  E=5.D-10 
B(NC)=H(NC) 
C(NC)=H(imC) 
B(NC+1  )=0. 
C(NC*1)=0,. 
NP=NC-1 

20  D049J=1,1000 
D021I 1=1, MP 
I=NC-I1 

B(I  )=HU)+R*B(  1  +  1) 

21  CU)=B(  I)  *R*C(  1*1) 
IF(DABS(B(1 )/H(l ) )-E)50,50,24 

24  IF(C(2))23, 22,23 

22  R=R+1. 
G0T030 

23  R=R-0<1)/C(  2) 

30  D037I 1=1, NP 
I =NC-I1 

B(  I  )  =  H(  I)-P*B(  H-l)-Q*B(I  +  2) 
37  C(I )  =  B( I)-P*C (I  +  1)-Q*C<  1  +  2) 
IF(H( 2J )32,31  ,32 

31  IF(DABS(B(2 )/H( I) )-E)33,33,34 
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32 
33 
34 


35 

36 
49 

50 

51 

52 
53 
54 

70 

71 

72 
73 
74 


75 

80 

81 
82 


76 
77 

100 


IF(DA 
IF(DA 
CBAR  = 
D=C(  3 
IF(D) 
P  =  P-2 
Q=Q*( 
GGT04 
P=P+( 
Q  =  G+( 
CONTI 
E  =  E*1 
G0T02 
NC=NC 
V(NC) 
I  F  (  IR 
U(NC) 
G0T05 
U(NC) 
D054I 
H(  I)  = 
G0TC4 
NC-NC 
IFtIR 
QP  =  1. 
PP  =  P/ 
G0T07 
QP=Q 
PP=P/ 
F  =  (PP 
IF(F) 
U  ( NC  + 
U(NC) 
V(NO 
V(NC) 
GGT07 
IF(PP 
U(NC  + 
GO  TO 
U(NO 
CONTI 
V(NO 
U  ( NC  ) 
V(NC) 
D077I 
H(  I)  = 
G0TD4 
RETUk 
END 


PROOT   (cont'd) 

eS(B(2)/H(2)5-E)33,33,34 
BS(B(1 )/H(l J )-F)70,70,34 
C(2)-B(2) 
)**2-CBAR*C<4) 

36,35, 36 


G  + 

9 

B( 

-B 

NU 

0. 

0 

-1 

=  0 

EV 

=  1 

3 

=R 

=  1 

B( 


1..J 

2)*C(3)-B(1)*C(4) )/D 
<2)*CBAR+B(1)*C(3)  )/D 

E 


)51,52,5? 
./R 


,NC 
1  +  1) 


-2 

EV>71,72,7? 

/Q 

(0*2.0) 

3 

2.0 

)**2-QP 

74,75,75 

i)=-PP 

=  -PP 

1)=DSQRT(-F) 

=-V<NC+l) 

6 

)81,30,81 

1)=-DSQRT<F) 

82 
1)=-(PP/DABS(PP) )*(DABS( PP)+DSQRT(F) ) 
NUE 
1)=0. 

=  QP/U(NOl) 
=  0. 
=  1,MC 
B( 1+2) 

N 
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DET 

FUNCTION    DET(A,KC) 
C  THIS    FUNCTION    SUBP^DGRAM    FINDS    THE    DETERMINANT    OF    A    MATRIX 

C  USING    DIAGOMALIZATION    PROCEDURE 

DOUBLE    PRECISIGN    A( I  0 , 1 0 )  ,B ( 10 , 10 ) , TEMP, DET 

IREV    =0 

D01    1=1, KC 

D01    J=l,KC 

1  B( I, J)=A( I, J) 
D020I=1,KC 
K=I 

9    IF(B(K,I) ) 10,11,10 

11  K=K+1 
IF(K-KC)    9,9,51 

10    IF(I-K)     12,14,51 

12  D013M=1 ,KC 
TEKP=ril I ,M) 
B( I,M)=B(K,M) 

13  B(K,M)=TEMP 
IREV    =     IREV+i 

14  11=1+1 

IFUI .GT.KC)    GO    TO    20 
D017    M=II,KC 

18  IF(B(M, I) )     19,17,19 

19  TEMP    =6(M,I  }/B(I  ,1) 
D016N=I,KC 

16  B(M,N)=B(M,N)-B< I ,N)*TEMP 

17  CONTINUE 

20  CONTINUE 
DET=1  . 
D02    1=1, KC 

2  DET=DET*B< 1,1) 
DET=(-1. )**IREV*DET 
RETURN 

51    DtT=0 
RETURN 
END 
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MPY 

SUBROUTINE    MPY ( 8 TG, N , ANS, NS ) 

DOUBLE    PRECISION!    3  (  1  0  J  ,G  I  10  i  ,  W  (  10  ,  10 )  ,  ANS<  I  L  )  ,  Z  ED  (  10  ,  10 ,  10 ) 
FORr'S    A    SCALAR    POLYNOMIAL    FROM    A    MATRIX    POLYNOMIAL 
COMMON    ZED 
DO    1       1=1, N 
DO    1     J  =  1,N 
W( J,  I  )=0.0 
DO    1    K=1,N 

1  W(  J, I  )=W(  J,  U+ZEDU,  J,K)*B(K) 
DO    2     I=1,N 

ANS( I )=0.0 
DO    2    J=1,N 

2  ANS( I)=ANS( I)*W( J,I)*G( J) 
NN=N*1 
ANS(NN)=0.0 
NS=N-1 
RETURN 
END 
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CHREQ 

SUBROUTINE  CHREQ1 A , N , C , NRM ) 
C   THIS  SUBROUTINE  FINDS  THE  COEFFICIENTS  OF  THE  CHARACTERISTIC 
C      POLYNOMIAL  USING  THE  LEVERRIER  ALGO°ITHM 

DOUBLE  PRECISION  A ( I  0 , 10 ) , C ( 1 1 ) , ATE MP( 10, 10 )  , PROO ( 10 , 1 0) , 
1  ZED( 10,10,10) 

COMMON  ZED 

DATA  ATEMP/100*0.0/ 

1000  FORMAT  ( 1H0  ,  5X  ,  31HT4E  MATRIX  COEFFICIENTS  3F  THE    , 
*    33HNUMERAT0R  OF  THE  RESOLVENT  MATRIX   ) 

1001  FORMAT  (1H0,   5X,29HTHE  MATRIX  COEFFICIENT  OF  S**,I1/) 

1002  FOPMAT  (1PoE20.7) 

1003  FORMAT  < 1H0  ,45( 1H* ) ) 
CALL  CHREQA(A,N,C) 
DO  6  5  I  =  1,N 

65  ATEMP( I, I)=1.0 

70  00  80  1*1 tN 
DO  30  J=1,N 

80  ZED(N,I,J)=ATEMP( I, J) 

IF  (NRM.NE.O)  GO  TO  71 

WRITE  (6,1003) 

WRITE  (6,1000) 

M=N-1 

WRITE  (6,1001)  M 

DO  3  5  I=1,M 
35  WRITE  (6,1002)  ( ATEMP ( I , J ) , J=l , N) 

71  DO  40  1*1, N 
DO  40  J=l ,N 

40  ATEMP( I , J) =A( I, J) 

DO  10  1=1, N 

NNN=N-I 

IF  (  I  .EQ.l )  GO  TO  55 

IF  4NRM.NE.0)  GO  TO  60 

WRITE  (6,1001)  NNN 

DO  45  J=l ,N 
45  WRITE  (6,1002)  ( ATEMP( J ,K ) , K=l , N) 
60  NP=NNN+1 

DO  90  11=1, N 

DO  90  J=1,N 
90  ZED(NP, II , J)=ATEMP(  II, J) 

DO  15  J=1,N 

DO  15  K=1,N 

PROD( J,K)=0.0 

DO  15  L=1,N 
15  PROD( J,K)=PROD(J,<)  +  ( A( J ,L) *ATEMP(  L , K ) ) 

DO  13  J=l ,N 

DO  13  K=1,N 
13  ATEM?(J,K)=PROD( J,K) 
55  DO  10  J=1,N 

NZ=N-I+1 
10  ATEMP(J,J)=ATEMP( J,J)+C(NZ) 

RETURN 

END 
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CHREQA 

SUBROUTINE    CHRF3A ( A . N T C ) 

DOUBLE    PRECISION    C ( 1 1 i  ,3( 1 0, 1 0 ) , A( 1 0 ,  10 ) , D( 300 ) 
DIMENSION    J (11) 
NN=N+1 

DO    20    1=1, NN 
20    C(  I)=0.0 
C(NN)=1.0 
DO    14    M=1,N 
K  =  0 
L  =  l 

J(  1 ) =1 
GO   TO    2 

1  J(L)=J(L)+1 

2  IF(L-M)    3,5,50 

3  MM=M-1 

DO    4    I=L,MM 
11=1+1 

4  J(  II  )=J( I)  +  l 

5  DO    10    1=1 , M 
DO    10    KK=1,M 
NR=J( I) 
NC=J(KK) 

10    B(  I,KK)=A(NR,NC) 
K  =  K+1 

D(K)=DET(B,M) 
DO    6     I=1,M 
L=M-I*1 
IF< J(L)-(N-M+L) )     1,6,50 

6  CONTINUE 
M1=N-M+1 

DO    14    !=1,K 
14    C(Ml)=C(Mi)+D(I)*(-1.0)**M 

RETURN 
50    PRINT    2000 
2000    FORMAT    < 1H0  ,  5X , i 5HERR0R    IN    CHREQA) 

RETURN 

END 
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DIMCH 

SUBROUTINE    DI MCH ( A, NA , 8 ) 
DOUBLE    PRECISION    A( 10  ,  10 ) , B ( 10, 10) 
DIMENSION    NA<2) 
II=NA( 1) 
JJ=NA(2) 
DO    1     1=1,11 
DO    1     J=1,JJ 
JARG=i J-1)*NA< 1 )*I 
6(1 ,J)=A( JARG,1) 
1       CONTINUE 
RETURN 
END 


A  35 


DISTRIBUTION  LIST 


No.  of  Copies 


Library 

Code  0212 

Naval  Postgraduate  School 

Monterey,  CA  939*4-0 


2.      Defense  Documentation  Center 
Cameron  Station 
Alexandria,  VA     223lU 


3.      Department  of  Aeronautics 
Code   57 

Naval  Postgraduate  School 
Prof.   R.  W.   Bell 
Asst.    Prof.   R.  A.   Hess 
Monterey,   CA     939^0 


1 
20 


h.      Dean  of  Research 
Code  023 

Naval  Postgraduate  School 
Monterey,  CA  939^0 


5.  Guidance  and  Navigation  Branch 
Flight  Systems  Research  Division 
Ames  Research  Center 

Mr.  George  Xenakis 

Mr.  Leonard  McGee 
Moffett  Field,  CA  9*4035 


1 
1 


U166501 


DUDLEY  KNOX  LIBRARY  -  RESEARCH  REPORTS 


5  6853  01057725  7 


